# Introduction
Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des operation statistique avancée grâce au langage R, à partir d’un jeu de données. Ainsi, notre choix s’est porté sur des données statistiques relatives à
Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer des régression linéaire sur notre jeu de données et produire un rapport.
# Outils et environnement de travail
R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.
La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.
L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.
Chargement des jeux de données :
data <- read_csv("data/prestige.csv") %>% as_tibble()
data
## # A tibble: 102 x 6
## education income women prestige census type
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 13.1 12351 11.2 68.8 1113 prof
## 2 12.3 25879 4.02 69.1 1130 prof
## 3 12.8 9271 15.7 63.4 1171 prof
## 4 11.4 8865 9.11 56.8 1175 prof
## 5 14.6 8403 11.7 73.5 2111 prof
## 6 15.6 11030 5.13 77.6 2113 prof
## 7 15.1 8258 25.6 72.6 2133 prof
## 8 15.4 14163 2.69 78.1 2141 prof
## 9 14.5 11377 1.03 73.1 2143 prof
## 10 14.6 11023 0.94 68.8 2153 prof
## # … with 92 more rows
En statistiques, et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.
Ici, nous testons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable éducation (qui reflète le niveau d’étude).
education <- data %>%
ggplot(aes(x = education, y = prestige)) +
geom_point() +
geom_smooth(method = loess, color = "red", fill = "#69b3a2", se = T) +
geom_smooth(method = lm, color = "blue", se = F) +
theme_ipsum()
scatter_edu <- ggplotly(education)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
scatter_edu
Nous testons aussi la relation entre la variable prestige et la variable income.
income <- data %>%
ggplot(aes(x = income, y = prestige)) +
geom_point() +
geom_smooth(method = loess, color = "red", fill = "#69b3a2", se = T) +
geom_smooth(method = lm, color = "blue", se = F) +
theme_ipsum()
scatter_inc <- ggplotly(income)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
scatter_inc
prest.lm1 <- lm(prestige~education, data = Prestige)
prest.lm1
##
## Call:
## lm(formula = prestige ~ education, data = Prestige)
##
## Coefficients:
## (Intercept) education
## -10.732 5.361
bacf <- acf(residuals(prest.lm1), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(prest.lm1)))
lag_plot <- bacfdf %>%
ggplot(aes(x = lag, y = acf)) +
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0)) +
geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +
geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
theme_ipsum()
lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot
durbinWatsonTest(prest.lm1)
## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0.006
## Alternative hypothesis: rho != 0
qq_edu <- prest.lm1 %>%
ggplot(aes(sample = residuals(prest.lm1) / 10))+
stat_qq() + stat_qq_line() +
#theme_ipsum() +
labs(x = "Theoretical Quantiles\nlm(prestige ~ education)",
y = "Standardized residuals")
#qq_edu <- ggplotly(qq_edu)
qq_edu
shapiro.test(residuals(prest.lm1))
##
## Shapiro-Wilk normality test
##
## data: residuals(prest.lm1)
## W = 0.98065, p-value = 0.1406
prest.lm2 <- lm(prestige~income, data = Prestige)
prest.lm2
##
## Call:
## lm(formula = prestige ~ income, data = Prestige)
##
## Coefficients:
## (Intercept) income
## 27.141176 0.002897
qq_inc <- prest.lm2 %>%
ggplot(aes(sample = residuals(prest.lm2) / 10))+
stat_qq() + stat_qq_line() +
#theme_ipsum() +
labs(x = "Theoretical Quantiles\nlm(prestige ~ income)",
y = "Standardized residuals")
#qq_inc <- ggplotly(qq_inc)
qq_inc
shapiro.test(residuals(prest.lm2))
##
## Shapiro-Wilk normality test
##
## data: residuals(prest.lm2)
## W = 0.97169, p-value = 0.02729
plot(prest.lm1, 3)
ncvTest(prest.lm1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545, Df = 1, p = 0.42635
plot(prest.lm1,1)
summary(prest.lm1)
##
## Call:
## lm(formula = prestige ~ education, data = Prestige)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16
confint(prest.lm1)
## 2.5 % 97.5 %
## (Intercept) -18.027220 -3.436744
## education 4.702223 6.019533
influenceIndexPlot(prest.lm1)
outlierTest(prest.lm1)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## newsboys -2.98896 0.0035306 0.36012
prest.lm1bis <- lm(prestige~education, data=Prestige[-c(53,67),])
compareCoefs(prest.lm1 ,prest.lm1bis)
## Calls:
## 1: lm(formula = prestige ~ education, data = Prestige)
## 2: lm(formula = prestige ~ education, data = Prestige[-c(53, 67), ])
##
## Model 1 Model 2
## (Intercept) -10.73 -11.26
## SE 3.68 3.53
##
## education 5.361 5.417
## SE 0.332 0.318
##
my_df <- data.frame(education=c(10.25))
predict(prest.lm1, newdata=my_df)
## 1
## 44.21701
predict(prest.lm1, newdata=my_df, interval="prediction")
## fit lwr upr
## 1 44.21701 26.06518 62.36885
predict(prest.lm1, newdata=my_df, interval="confidence")
## fit lwr upr
## 1 44.21701 42.40008 46.03395
my_df <- data.frame(education=c(10.25, 11.25, 12.25))
predict(prest.lm1, newdata=my_df,interval="confidence")
## fit lwr upr
## 1 44.21701 42.40008 46.03395
## 2 49.57789 47.75810 51.39768
## 3 54.93877 52.89190 56.98564
my_pres <- Prestige
my_pres$res <-residuals(prest.lm1)
head(my_pres)
## education income women prestige census type res
## gov.administrators 13.11 12351 11.16 68.8 1113 prof 9.250875
## general.managers 12.26 25879 4.02 69.1 1130 prof 14.107621
## accountants 12.77 9271 15.70 63.4 1171 prof 5.673573
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof 6.310758
## chemists 14.62 8403 11.68 73.5 2111 prof 5.855950
## physicists 15.64 11030 5.13 77.6 2113 prof 4.487854
my_pres$fitted <-fitted(prest.lm1)
head(my_pres)
## education income women prestige census type res
## gov.administrators 13.11 12351 11.16 68.8 1113 prof 9.250875
## general.managers 12.26 25879 4.02 69.1 1130 prof 14.107621
## accountants 12.77 9271 15.70 63.4 1171 prof 5.673573
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof 6.310758
## chemists 14.62 8403 11.68 73.5 2111 prof 5.855950
## physicists 15.64 11030 5.13 77.6 2113 prof 4.487854
## fitted
## gov.administrators 59.54913
## general.managers 54.99238
## accountants 57.72643
## purchasing.officers 50.48924
## chemists 67.64405
## physicists 73.11215
head(predict(prest.lm1))
## gov.administrators general.managers accountants purchasing.officers
## 59.54913 54.99238 57.72643 50.48924
## chemists physicists
## 67.64405 73.11215
head(fitted(prest.lm1))
## gov.administrators general.managers accountants purchasing.officers
## 59.54913 54.99238 57.72643 50.48924
## chemists physicists
## 67.64405 73.11215
vcov(prest.lm1)
## (Intercept) education
## (Intercept) 13.520978 -1.1835055
## education -1.183505 0.1102162
ggplot(Prestige, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)")
## `geom_smooth()` using formula 'y ~ x'
int_pred <- predict(prest.lm1, interval="prediction")
## Warning in predict.lm(prest.lm1, interval = "prediction"): predictions on current data refer to _future_ responses
my_prest2 <-cbind(Prestige, int_pred)
head(my_prest2)
## education income women prestige census type fit
## gov.administrators 13.11 12351 11.16 68.8 1113 prof 59.54913
## general.managers 12.26 25879 4.02 69.1 1130 prof 54.99238
## accountants 12.77 9271 15.70 63.4 1171 prof 57.72643
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof 50.48924
## chemists 14.62 8403 11.68 73.5 2111 prof 67.64405
## physicists 15.64 11030 5.13 77.6 2113 prof 73.11215
## lwr upr
## gov.administrators 41.33302 77.76523
## general.managers 36.81573 73.16903
## accountants 39.52816 75.92469
## purchasing.officers 32.33470 68.64379
## chemists 49.31584 85.97226
## physicists 54.67820 91.54609
ggplot(my_prest2, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
geom_line(aes(y=lwr), color = "red", linetype = "dashed")+
geom_line(aes(y=upr), color = "red", linetype = "dashed")+
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)")
## `geom_smooth()` using formula 'y ~ x'